* Hafstead-Williams-Chen Search-CGE model

* Produces results for central calibration carbon tax with lump-sum rebates

* Define sets 

sets
        t "time periods" /t1*t51/
        tp(t) /t1*t51/
        tpi(t) /t1/
        tfirst(t) "First time period"
        tlast(t) "Final time period"
        ind /i1*i22/
        ie(ind) "Energy Industries" /i1*i2,i5*i7/
        im(ind) "Material Industries" /i3,i4,i8*i22/
        oilgas(ind) "Oil&Gas Extraction" /i1/
        coal(ind) "Coal Mining" /i2/
        industrial(ind) "Industrial" /i1*i4,i8*i12/
        commercial(ind) "Commercial" /i13*i14,i16*i21/
        gov(ind) "Government" /i22/
        transportation(ind) "Transportation" /i15/
        elec(ind) "Elecric Power" /i5/ 
        ngd(ind) "Nat Gas Dist" /i6/
        refine(ind) "Refine" /i7/
        refine_n(ind) /i1*i6,i8*i22/
        primary(ind) "Primary Energy" /i1*i2/ 
        secondary(ind) "Secondary Energy" /i6*i7/ 
        durable(ind) "Durable Manuf" /i11/ 
        nondurable(ind)  "Non-durable Manuf" /i12/
        other(ind) /i3,i4,i8,i9,i10,i16/
        trade(ind) "Trade"  /i13,i14/
        it "iterations" /it1*it8/
        out /out1*out12/
        services(ind) /i17*i21/
        i /i1*i22,Gov,All/  
;

tfirst(t) = yes$(ord(t) eq 1);
tlast(t) = yes$(ord(t) eq card(t));
alias (ind,ind2)
alias (ind,ind3)
alias(ie,ie2);
alias(im,im2);
alias(industrial,industrial2)
alias(commercial,commercial2)
alias(refine_n,refine_n2)

scalar
nt "Number of Periods" /43/
;


* DEFINE DEEP PARAMETERS
parameters
         beta "time preference" /.996/
         betat(t)
         chi "Frisch Elasticity of Labor Supply" /1/
         cost_fire "Firing Costs" /.001/
         delta(*,*)
         tau "Nash bargaining weight" /.5/
         pi(ind) "Exogenous quit rate"
         pig "Government exogenous quit rate"
         pi_f(ind) "Foreign exogenous quit rate"
         pig_f "Foreign government exogenous quit rate"
         gamma "Curvature of matching function" /0.5/
         mu_e(ind,ind2) "Emissions factor"
         tau_l0 "Tax on Labor Income" /0.31/
         tau_p0 "Payroll tax" /0.06/
         sigma /1/
         sigmac "Elast. of subs. between cons. goods" /0.75/
         sig_e(ind) "Elast. of subs. between energy goods"
         sig_eg "Elast. of subs. between energy goods - gov"
         sig_m(ind) "Elast. of subs. between non-energy goods"
         sig_mg "Elast. of subs. between non-energy goods - gov"
         sig_v(ind) "Elast. of subs. between energy and non-energy goods"
         sig_vg "Elast. of subs. between energy and non-energy goods - gov"
         sig_y(ind) "Elast. of subs. between labor and inter. inputs"
         sig_yg "Elast. of subs. between labor and inter. inputs - gov"
         rho(ind)
         rhog
        rhos
        rhosg


         sig_df(ind)
         sig_c(ind)
         sig_c_f(ind)

;

* Sensitivity: increase match elasticity to 0.75

gamma = 0.75;

betat(t) = beta**(ord(t)-1);

sig_y(ind) = 0.5;
sig_e(ind) = 0.4;
sig_m(ind) = 0.6;
sig_v(ind) = 0.7;

sig_v('i5')  = 0.4;
sig_v('i6')  = 0.4;
sig_v('i7')  = 0.4;
sig_v('i8')  = 0.4;

sig_eg = 0.2;
sig_mg = 0.2;
sig_vg = 0.2;
sig_yg = 0.2;

rho(ind) = (sig_y(ind)-1)/sig_y(ind);
rhog = (sig_yg-1)/sig_yg;

sig_e(elec) = 1.5;

* Foreign and domestic sustituion elas for both production and consumption
sig_df(ind) = 1.5;
sig_c(ind) = 1.5;
sig_c_f(ind) = 1.5;

rhos(ind) = 0;
rhosg = 0;

delta(ind,ind2) = 0;
delta(ind,ind2)$(ord(ind2) = ord(ind)) = 1;



$INCLUDE HWC_JAERE_CAL_QFIX.gms



* DEFINE STEADY STATE CARBON PRICE
parameter
p_e0
tar_i2(ind)
gammay_ss_ct(ind)
gammay(ind,t)
gammav_ss(ind)
gammav(ind,t)
;

p_e0 = 0;
tar_i2(refine) = f_y0(refine)*sum(oilgas,mu_e(oilgas,refine));
gammay_ss_ct(ind) = gammay0.l(ind);
gammay(ind,t) = gammay0.l(ind);
gammav_ss(ind) = gammav0(ind);
gammav(ind,t) = gammav0(ind);


$INCLUDE HWC_JAERE_EQN_SS.gms



options decimals = 5;




* REPLICATE BASE CASE WITH P_E = 0;


options iterlim = 0;

p_e_ss(ind,ind2) = 0;
p_e_c_ss(ind) = 0;
sub_ss_ct.fx(ind) = 0;


delta_tau.fx = 1;


solve hwc_ss_ct using mcp ;



* SOLVE STEADY STATE WITH P_E>0
options iterlim = 100000000;

* Define steady state output parameters

parameter
utot_ss_ct
n_d
n_d_perc
n_ss_it
n_ss_ct_it

l_ss_ct
l_d
l_d_perc
h_d
h_d_perc
ng_d
ng_d_perc
lg_ss_ct
lg_d
lg_d_perc
hg_d
hg_d_perc
v_ss
v_ss_ct
v_d
v_d_perc
vg_ss
vg_ss_ct
vg_d
vg_d_perc
tlump_d
net_gross
crev_ss_ct
perc_e
perc_e_elec
perc_e_y
perc_l
perc_h
perc_n
perc_v

tbal_ss_db
mc_ss_db
bc_ss_db
mc_ss_f_db
bc_ss_f_db
yg_ss_db
yg_ss_f_db
sub_ss_db

tbal_db
mc_db
bc_db
mc_f_db
bc_f_db
yg_db
yg_f_db

ndiff
ngdiff
nfdiff
nfgdiff

utot
utot_f
utot_diff
sub_it
tau_l_d

perc_l_new
perc_h_new
perc_n_new

vbar_perc
vbarg_perc
perc_earn
earn_perc
w_perc
perc_w
perc_lp
perc_income

p_e_it

hours_db1
hours_db2
hours_db3
hours_db4
hours_db5

b_rr_ss_ct
p_f_ss_perc
;

* Define transition output parameters

p_e0 = 0;
loop(it,p_e0 = p_e0+5;

p_e_ss(ind,ind2) = p_e0;
p_e_c_ss(ind) = p_e0;
p_e_it(it) = p_e0;



solve hwc_ss_ct using mcp   ;




* SAVE STEADY STATE OUTPUT

b_rr_ss_ct(it) = (b0.l)/(( (sum(ind,(1-tau_l_ss_ct.l)*w_ss_ct.l(ind)*h_ss_ct.l(ind)*n_ss_ct.l(ind)) + (1-tau_l_ss_ct.l)*wg_ss_ct.l*hg_ss_ct.l*ng_ss_ct.l)/(sum(ind,n_ss_ct.l(ind)) + ng_ss_ct.l)  )/p_cbar_ss_ct.l);

utot_ss_ct(it) = 1-ntot_ss_ct.l;

n_d(ind,it)= n_ss_ct.l(ind)-n_ss.l(ind);
n_d_perc(ind,it)= (n_ss_ct.l(ind)-n_ss.l(ind))/n_ss.l(ind);

n_ss_it(ind,it) = n_ss.l(ind);
n_ss_ct_it(ind,it) = n_ss_ct.l(ind);

ng_d(it)= ng_ss_ct.l-ng_ss.l;
ng_d_perc(it)= (ng_ss_ct.l-ng_ss.l)/ng_ss.l;

l_ss_ct(ind) = h_ss_ct.l(ind)*n_ss_ct.l(ind);
*(1-vbar_ss_ct.l(ind));
l_d(ind) = l_ss_ct(ind)-l_ss(ind);
l_d_perc(ind,it) = (l_ss_ct(ind)-l_ss(ind))/l_ss(ind);

lg_ss_ct = hg_ss_ct.l*ng_ss_ct.l;
*(1-vbarg_ss_ct.l);
lg_d = lg_ss_ct-lg_ss;
lg_d_perc(it) = (lg_ss_ct-lg_ss)/lg_ss;

h_d(ind,it) = h_ss_ct.l(ind) - h_ss.l(ind);
h_d_perc(ind,it) = (h_ss_ct.l(ind) - h_ss.l(ind))/h_ss.l(ind);

hg_d(it) = hg_ss_ct.l - hg_ss.l;
hg_d_perc(it) = (hg_ss_ct.l - hg_ss.l)/hg_ss.l;

v_ss(ind) = h_ss.l(ind)*n_ss.l(ind)*vbar_ss.l(ind);
v_ss_ct(ind) = h_ss_ct.l(ind)*n_ss_ct.l(ind)*vbar_ss_ct.l(ind);
v_d(ind) = v_ss_ct(ind)-v_ss(ind);
v_d_perc(ind,it) = v_d(ind)/v_ss(ind);

vg_ss = hg_ss.l*ng_ss.l*vbarg_ss.l;
vg_ss_ct = hg_ss_ct.l*ng_ss_ct.l*vbarg_ss_ct.l;
vg_d = vg_ss_ct-vg_ss;
vg_d_perc(it) = vg_d/vg_ss;

vbar_perc(ind,it) = vbar_ss_ct.l(ind)/vbar_ss.l(ind) - 1;
vbarg_perc(ind,it) = vbarg_ss_ct.l/vbarg_ss.l - 1;


tlump_d(it) = tlump_ss_ct.l - tlump_ss;

earn_perc(ind,it) = (h_ss_ct.l(ind)*w_ss_ct.l(ind)/p_cbar_ss_ct.l)/(h_ss.l(ind)*w_ss.l(ind)/p_cbar_ss) - 1;
w_perc(ind,it) = (w_ss_ct.l(ind)/p_cbar_ss_ct.l)/(w_ss.l(ind)/p_cbar_ss) - 1;


perc_e(it) = (etot_ss_ct.l-etot0)/etot0;
perc_e_elec(elec,it) = (e_ss_ct.l(elec) - e0(elec))/e0(elec);

perc_e_y(elec,it) = (e_ss_ct.l(elec)/y_ss_ct.l(elec) - e_y0(elec))/e_y0(elec);
perc_l(it) = (sum(ind,l_ss_ct(ind))+lg_ss_ct-sum(ind,l_ss(ind))-lg_ss)/(sum(ind,l_ss(ind))+lg_ss);
perc_h(it) = (sum(ind,n_ss.l(ind)*h_ss_ct.l(ind))+ng_ss.l*hg_ss_ct.l-sum(ind,n_ss.l(ind)*h_ss.l(ind))-ng_ss.l*hg_ss.l)/(sum(ind,n_ss.l(ind)*h_ss.l(ind))+ng_ss.l*hg_ss.l);
perc_n(it) = (sum(ind,n_ss_ct.l(ind))+ng_ss_ct.l-sum(ind,n_ss.l(ind))-ng_ss.l)/(sum(ind,n_ss.l(ind))+ng_ss.l);
perc_v(it) = (sum(ind,v_ss_ct(ind))+vg_ss_ct-sum(ind,v_ss(ind))-vg_ss)/(sum(ind,v_ss(ind))+vg_ss);
perc_lp(it) = (sum(ind,n_ss_ct.l(ind)*h_ss_ct.l(ind)*(1-vbar_ss_ct.l(ind))) + ng_ss_ct.l*hg_ss_ct.l*(1-vbarg_ss_ct.l))/(sum(ind,n_ss.l(ind)*h_ss.l(ind)*(1-vbar_ss.l(ind))) + ng_ss.l*hg_ss.l*(1-vbarg_ss.l)) - 1   ;
perc_w(it) = (  (sum(ind,n_ss.l(ind)*w_ss_ct.l(ind))+ng_ss.l*wg_ss_ct.l)/p_cbar_ss_ct.l-(sum(ind,n_ss.l(ind)*w_ss.l(ind))+ng_ss.l*wg_ss.l)/p_cbar_ss )/( (sum(ind,n_ss.l(ind)*w_ss.l(ind))+ng_ss.l*wg_ss.l)/p_cbar_ss);
perc_earn(it) = (  (sum(ind,n_ss.l(ind)*w_ss_ct.l(ind)*h_ss_ct.l(ind))+ng_ss.l*wg_ss_ct.l*hg_ss_ct.l)/p_cbar_ss_ct.l-(sum(ind,n_ss.l(ind)*w_ss.l(ind)*h_ss.l(ind))+ng_ss.l*wg_ss.l*hg_ss.l)/p_cbar_ss    )   / ( (sum(ind,n_ss.l(ind)*w_ss.l(ind)*h_ss.l(ind))+ng_ss.l*wg_ss.l*hg_ss.l)/p_cbar_ss);
perc_income(it) = (  (sum(ind,n_ss.l(ind)*w_ss_ct.l(ind)*h_ss_ct.l(ind))+ng_ss.l*wg_ss_ct.l*hg_ss_ct.l + (1-ntot_ss_ct.l)*p_cbar_ss_ct.l*b0.l)/p_cbar_ss_ct.l-(sum(ind,n_ss.l(ind)*w_ss.l(ind)*h_ss.l(ind))+ng_ss.l*wg_ss.l*hg_ss.l + (1-ntot_ss.l)*p_cbar_ss*b0.l)/p_cbar_ss    )   / ( (sum(ind,n_ss.l(ind)*w_ss.l(ind)*h_ss.l(ind))+ng_ss.l*wg_ss.l*hg_ss.l + (1-ntot_ss.l)*p_cbar_ss*b0.l)/p_cbar_ss);



tbal_ss_db(it) = sum(ind,p_ss_ct.l(ind)*exports_ss_ct.l(ind) - (p_f_ss_ct.l(ind)/exch_ss_ct.l)*imports_ss_ct.l(ind)) ;
mc_ss_db(ind,it) = y_ss_ct.l(ind) - cd_ss_ct.l(ind) - gd_ss_ct.l(ind) -  exports_ss_ct.l(ind) - sum(ind2,iod_ss_ct.l(ind,ind2)) ;
bc_ss_db(it) = (1-tau_l_ss_ct.l)*(sum(ind,n_ss_ct.l(ind)*w_ss_ct.l(ind)*h_ss_ct.l(ind)) + ng_ss_ct.l*wg_ss_ct.l*hg_ss_ct.l) + (1-ntot_ss_ct.l)*p_cbar_ss_ct.l*b0.l + p_cbar_ss_ct.l*tlump_ss_ct.l  + sum(ind,prof_ss_ct.l(ind))  - p_cbar_ss_ct.l*cbar_ss_ct.l;
mc_ss_f_db(ind,it) = y_f_ss_ct.l(ind) - cd_f_ss_ct.l(ind) - gd_f_ss_ct.l(ind) -  imports_ss_ct.l(ind) - sum(ind2,iod_f_ss_ct.l(ind,ind2)) ;
bc_ss_f_db(it) = (1-tau_l0)*(sum(ind,n_f_ss_ct.l(ind)*w_f_ss_ct.l(ind)*h_f_ss_ct.l(ind)) + ng_f_ss_ct.l*wg_f_ss_ct.l*hg_f_ss_ct.l) + (ratiof-ntot_f_ss_ct.l)*p_cbar_f_ss_ct.l*b0_f.l + p_cbar_f_ss_ct.l*tlump_f_ss_ct.l  + sum(ind,prof_f_ss_ct.l(ind))  - p_cbar_f_ss_ct.l*cbar_f_ss_ct.l;
yg_ss_db(it) =  gammayg0.l*( alphapg.l*(hg_ss_ct.l*ng_ss_ct.l*(1-vbarg_ss_ct.l))**rhog + (1-alphapg.l)*vvg_ss_ct.l**rhog)**(1/rhog) - yg0;
yg_ss_f_db(it) =  gammayg0_f.l*( alphapg_f.l*(hg_f_ss_ct.l*ng_f_ss_ct.l*(1-vbarg_f_ss_ct.l))**rhog + (1-alphapg_f.l)*vvg_f_ss_ct.l**rhog)**(1/rhog) - yg0_f;
sub_ss_db(ind,it) = sum(ie,p_e_ss(ie,ind)*mu_e(ie,ind)*(iod_ss_ct.l(ie,ind)+iof_ss_ct.l(ie,ind)) + tar_i2(ie)*p_e_ss(ie,ind)*iof_ss_ct.l(ie,ind) )- p_ss_ct.l(ind)*sub_ss_ct.l(ind)*y_ss_ct.l(ind);

p_f_ss_perc(ind,it) = p_f_ss_ct.l(ind)/exch_ss_ct.l;

if(ord(it) = 8,
hours_db1(ind,it) = chi*(log(h_ss_ct.l(ind))-log(h_ss.l(ind)));
hours_db2(ind,it) = (log(lambda_ss_ct.l)-log(lambda_ss)) + (log(f_l_ss_ct.l(ind))-log(f_l_ss.l(ind))) + (log(p_ss_ct.l(ind))-log(p_ss(ind)));
hours_db3(ind,it) = (log(lambda_ss_ct.l)-log(lambda_ss));
hours_db4(ind,it) = (log(f_l_ss_ct.l(ind))-log(f_l_ss.l(ind))) + (log(p_ss_ct.l(ind))-log(p_ss(ind))); ;
hours_db5(ind,it) = (log(p_ss_ct.l(ind))-log(p_ss(ind)));


);

);


* DISPLAY KEY RESULTS

display
utot_ss_ct
n_d_perc
l_d_perc
h_d_perc
ng_d_perc
lg_d_perc
hg_d_perc
perc_e_y
perc_e_elec
perc_e
perc_n
perc_h
perc_l
vbar_perc
vbarg_perc
tau_l_ss_ct.l
earn_perc
w_perc


tbal_ss_db
mc_ss_db
bc_ss_db
mc_ss_f_db
bc_ss_f_db
yg_ss_db
yg_ss_f_db
sub_ss_db

hours_db1
hours_db2
hours_db3
hours_db4
hours_db5

b_rr_ss_ct
p_f_ss_perc



;

parameter
table_HWC_LUMP_p_e40_ME1
;


table_HWC_LUMP_p_e40_ME1('out1') = utot_ss_ct('it8')-.05;


execute_unload "table_JAERE_HWC_LUMP_p_e40_ME1.gdx"     table_HWC_LUMP_p_e40_ME1

display

table_HWC_LUMP_p_e40_ME1
;

